--****
-- from ck in  boggle.zip http://www.usingeuphoria.com
--
-- NOTES:
-- Some routines use recursion, and may be slow and memory intensive for large sequences
--
-- A more sophisticated method is needed if permutations or selections
-- involving sets longer than 13 are needed.
-- The sequence factorials is only 13 long because 13! is the biggest
-- that can be stored as a single element without loss of accuracy.
--
-- Permutations and Combinations of up to 13 things
-- #  permutation q of s: - perm_qs(3,"abcd") gives "acbd"
-- # selection q of r things from s: - select_qrs(19,3,"abcde") gives "bde"
-- # permutation q of r things from s: - perm_qrs(15,3,"abcd") gives "cad"
-- 
-- the ordering method utilized allows us to calculate 
-- a permutation or combination directly given its number (q).
-- * "perm_qs()"  returns permutation number q of the sequence s.
-- * "select_qrs()" returns combination number q of r elements from the sequence s.
-- * "perm_qrs()" returns permutation number q of r elements from the sequence s.
-- * There are no restrictions in the type or length in the elements of s, 
--    but, the length of s must not be more than 13.
-- 
-- changes ne1 11/09\\
-- :augmented with info/code from  
-- pynms Keith Dart <kdart@kdart.com>
-- License: LGPL  combinatorics.py,v 1.1 2004/01/09 02:52:28
-- slightly different nCr and added nPr and get_*\\
-- nCr is also similar in Don Cahela's math.e \\
-- maybe add Cartesian like in probstat which wraps gsl for python 
--
-- using bigfixedmath from Unkmar for calculations on n>13
--
-- permutions and combinations with repitions not currently available?
--testing is ongoing, 
--routines may well return less than all possible permutions and combinations

namespace permcomb

--****
-- The file "examples/trypc.ex" shows some use of "permcomb.e".
--
-- <eucode>
--  sequence word = "", letters="ABC"
-- 	for q=1 to nPr(3,2) do
-- 		word &= perm_qrs( q, 2, letters )
-- 	end for
--  printf(1, "%s\n", {word})
-- --ABBAACCAAB
-- </eucode>
-- 

--~ include std/math.e
public include std/sequence.e
public include qt/bigfixedmath.e
include std/text.e

--**
--definitive combinations from gsl  1..n 2..k
-- can use as indexes to your sequence <=length(6) , base 0
export constant combsu = {
 {
   {  1,  1,{0}}
 }
 ,{
    {  2,  2,{{0,1}}}
 }
 ,{
   {3,  2,{{0,1},{0,2},{1,2}}}
 , { 3,  3,{{0,1,2}}}
 }
 ,{{  4,  2,{{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}}
 ,{  4,  3,{{0,1,2},{0,1,3},{0,2,3},{1,2,3}}}
 ,{  4,  4,{{0,1,2,3}}}
 }
 ,{{  5,  2,{{0,1},{0,2},{0,3},{0,4},{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}}}
 ,{  5,  3,{{0,1,2},{0,1,3},{0,1,4},{0,2,3},{0,2,4},{0,3,4},{1,2,3},{1,2,4},{1,3,4},{2,3,4}}}
 ,{  5,  4,{{0,1,2,3},{0,1,2,4},{0,1,3,4},{0,2,3,4},{1,2,3,4}}}
 ,{  5,  5,{{0,1,2,3,4}}}
 }
 ,{{  6,  2,{{0,1},{0,2},{0,3},{0,4},{0,5},{1,2},{1,3},{1,4},{1,5},{2,3},{2,4},{2,5},{3,4},{3,5},{4,5}}}
 ,{  6,  3,{{0,1,2},{0,1,3},{0,1,4},{0,1,5},{0,2,3},{0,2,4},{0,2,5},{0,3,4},{0,3,5},{0,4,5},{1,2,3},{1,2,4},{1,2,5},{1,3,4},{1,3,5},{1,4,5},{2,3,4},{2,3,5},{2,4,5},{3,4,5}}}
 ,{  6,  4,{{0,1,2,3},{0,1,2,4},{0,1,2,5},{0,1,3,4},{0,1,3,5},{0,1,4,5},{0,2,3,4},{0,2,3,5},{0,2,4,5},{0,3,4,5},{1,2,3,4},{1,2,3,5},{1,2,4,5},{1,3,4,5},{2,3,4,5}}}
 ,{  6,  5,{{0,1,2,3,4},{0,1,2,3,5},{0,1,2,4,5},{0,1,3,4,5},{0,2,3,4,5},{1,2,3,4,5}}}
 ,{  6,  6,{{0,1,2,3,4,5}}}
 } }


--**
-- All permutations of 1..5 from gsl
-- can use as indexes to your {{{sequence, <length(6) base 0}}}
export constant permsu = {       
  {  1,{{0}}
  }
 ,{ 2,{{0,1},{1,0}}
  }
 ,{  3,{{0,1,2},{0,2,1},{1,0,2},{1,2,0},{2,0,1},{2,1,0}}
  }--1/4k
 ,{  4, {{0,1,2,3},{0,1,3,2},{0,2,1,3},{0,2,3,1},{0,3,1,2},{0,3,2,1},{1,0,2,3},{1,0,3,2},{1,2,0,3},{1,2,3,0},{1,3,0,2},{1,3,2,0},{2,0,1,3},{2,0,3,1},{2,1,0,3},{2,1,3,0},{2,3,0,1},{2,3,1,0},{3,0,1,2},{3,0,2,1},{3,1,0,2},{3,1,2,0},{3,2,0,1},{3,2,1,0}}
  }--over 1k
 ,{  5, {{0,1,2,3,4},{0,1,2,4,3},{0,1,3,2,4},{0,1,3,4,2},{0,1,4,2,3},{0,1,4,3,2},{0,2,1,3,4},{0,2,1,4,3},{0,2,3,1,4},{0,2,3,4,1},{0,2,4,1,3},{0,2,4,3,1},{0,3,1,2,4},{0,3,1,4,2},{0,3,2,1,4},{0,3,2,4,1},{0,3,4,1,2},{0,3,4,2,1},{0,4,1,2,3},{0,4,1,3,2},{0,4,2,1,3},{0,4,2,3,1},{0,4,3,1,2},{0,4,3,2,1},{1,0,2,3,4},{1,0,2,4,3},{1,0,3,2,4},{1,0,3,4,2},{1,0,4,2,3},{1,0,4,3,2},{1,2,0,3,4},{1,2,0,4,3},{1,2,3,0,4},{1,2,3,4,0},{1,2,4,0,3},{1,2,4,3,0},{1,3,0,2,4},{1,3,0,4,2},{1,3,2,0,4},{1,3,2,4,0},{1,3,4,0,2},{1,3,4,2,0},{1,4,0,2,3},{1,4,0,3,2},{1,4,2,0,3},{1,4,2,3,0},{1,4,3,0,2},{1,4,3,2,0},{2,0,1,3,4},{2,0,1,4,3},{2,0,3,1,4},{2,0,3,4,1},{2,0,4,1,3},{2,0,4,3,1},{2,1,0,3,4},{2,1,0,4,3},{2,1,3,0,4},{2,1,3,4,0},{2,1,4,0,3},{2,1,4,3,0},{2,3,0,1,4},{2,3,0,4,1},{2,3,1,0,4},{2,3,1,4,0},{2,3,4,0,1},{2,3,4,1,0},{2,4,0,1,3},{2,4,0,3,1},{2,4,1,0,3},{2,4,1,3,0},{2,4,3,0,1},{2,4,3,1,0},{3,0,1,2,4},{3,0,1,4,2},{3,0,2,1,4},{3,0,2,4,1},{3,0,4,1,2},{3,0,4,2,1},{3,1,0,2,4},{3,1,0,4,2},{3,1,2,0,4},{3,1,2,4,0},{3,1,4,0,2},{3,1,4,2,0},{3,2,0,1,4},{3,2,0,4,1},{3,2,1,0,4},{3,2,1,4,0},{3,2,4,0,1},{3,2,4,1,0},{3,4,0,1,2},{3,4,0,2,1},{3,4,1,0,2},{3,4,1,2,0},{3,4,2,0,1},{3,4,2,1,0},{4,0,1,2,3},{4,0,1,3,2},{4,0,2,1,3},{4,0,2,3,1},{4,0,3,1,2},{4,0,3,2,1},{4,1,0,2,3},{4,1,0,3,2},{4,1,2,0,3},{4,1,2,3,0},{4,1,3,0,2},{4,1,3,2,0},{4,2,0,1,3},{4,2,0,3,1},{4,2,1,0,3},{4,2,1,3,0},{4,2,3,0,1},{4,2,3,1,0},{4,3,0,1,2},{4,3,0,2,1},{4,3,1,0,2},{4,3,1,2,0},{4,3,2,0,1},{4,3,2,1,0}}
  }-- 6 is over 10k etc
}
--eventually will make a database with an outside program 
--accessabvle with more sane indexes and presorted
--same with permutations. if nothing else, then as test data
--though will save time on longer permutes 
--you get higher lenghts practically for free once the program works
-- init_combperm()

--**
--precalculated nPr sequence <21
export sequence nPrs={
 "0", "1", "4", "15", "64", "325", "1956", "13699", "109600", "986409"
 ,"9864100", "108505111", "1302061344", "16926797485"
 ,"236975164804", "3554627472075", "56874039553216"
 ,"966858672404689", "17403456103284420", "330665665962403999"
}

--**
--precalculated Fibonacci sequence <28
export sequence fibs={
--  0, 1
 {1, 1 }         ,{2, 1 }         ,{3, 2 }
 ,{4, 3 }        ,{5, 5 }         ,{6, 8 }
 ,{7, 13 }      ,{8, 21 }       ,{9, 34 }
 ,{10, 55 }    ,{11, 89 }     ,{12, 144 }
 ,{13, 233 }   ,{14, 377 }   ,{15, 610 }
 ,{16, 987 }   ,{17, 1597 }   ,{18, 2584 }
 ,{19, 4181 }  ,{20, 6765 }  ,{21, 10946 }
 ,{22, 17711 } ,{23, 28657} ,{24, 46368} 
 ,{25, 75025}  ,{26, 121393} ,{27, 196418}
 --~ ,{28, 317811} ,{29, 514229} ,{30, 832040}
 --~ ,{31, 1346269}
  --~ 2178309, 3524578, 5702887, 9227465,
--~ 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296,
--~ 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976,
--~ 7778742049, 12586269025, 20365011074, 32951280099, 53316291173,
--~ 86267571272, 139583862445, 225851433717, 365435296162, 591286729879,
--~ 956722026041, 1548008755920, 2504730781961, 4052739537881,
--~ 6557470319842, 10610209857723, 17167680177565, 27777890035288,
--~ 44945570212853, 72723460248141, 117669030460994, 190392490709135,
--~ 308061521170129, --71
}


--**
--precalculated factorials sequence <14
export constant factorials =  {
     1,2,6,24,120,720,5040,40320,362880,3628800
     ,39916800,479001600,6227020800
     }

sequence facts = { --~ Reach of computer numbers.
                                {"1", "1!"}
                                ,{"2", "2!"}
                                ,{"6", "3!"}
                               ,{"24", "4!"}
                              ,{"120", "5!   8-bit unsigned"}
                              ,{"720", "6!"}
                             ,{"5040", "7!"}
                            ,{"40320", "8!  16-bit unsigned"}
                           ,{"362880", "9!"}
                          ,{"3628800", "10!"}
                         ,{"39916800", "11!"}
                        ,{"479001600", "12!  32-bit unsigned"}
                       ,{"6227020800", "13!"}
                      ,{"87178291200", "14!"}
                    ,{"1307674368000", "15!"}
                   ,{"20922789888000", "16!"}
                  ,{"355687428096000", "17!"}
                 ,{"6402373705728000", "18!"}
               ,{"121645100408832000", "19!"}
              --~ ,{"2432902008176640000", "20!  64-bit unsigned"}
              ,{"2.43290200817664e+18", "20!  64-bit unsigned"}
             --~ ,{"51090942171709440000", "21!"}
             ,{"5.109094217170944e+19", "21!"}
           --~ ,{"1124000727777607680000", "22!"}
           ,{"1.12400072777760768e+21", "22!"}
          --~ ,{"25852016738884976640000", "23!"}
          ,{"25852016738884976640000", "23!"}
         --~ ,{"620448401733239439360000", "24!"}
         ,{"620448401733239439360000", "24!"}
       --~ ,{"15511210043330985984000000", "25!"}
       ,{"15511210043330985984000000", "25!"}
      --~ ,{"403291461126605635584000000", "26!"}
      ,{"403291461126605635584000000", "26!"}
    --~ ,{"10888869450418352160768000000", "27!"}
    ,{"10888869450418352160768000000", "27!"}
   --~ ,{"304888344611713860501504000000", "28!"}
   ,{"304888344611713860501504000000", "28!"}
  --~ ,{"8841761993739701954543616000000", "29!"}
  ,{"8841761993739701954543616000000", "29!"}
--~ ,{"265252859812191058636308480000000", "30!"}
--~ ,{"265252859812191058636308480000000", "30!"}
}


--~ bigprime

 --memoize bigfib results use append instead of swap
 --little problematic, would have to generate between the gaps too
--~ sequence prev{0,1} 

--**
-- non recursive bigfixed version.
--http://en.wikipedia.org/wiki/Fibonacci_number
--this would normally only be used if result could go over integer size
--returned a bigfixed number
export function bigfib(integer n)
    --make version with iadd instead
    --check with a a bigfixed too
  --~ bfi:set_places(1)   --is this local set 
  integer fibslen = length(fibs)
   --for <n30? return answers from a sequence
  if n <= 0 then
        return text_num(sprint(0)) 
  elsif n <= fibslen then
        return text_num(sprint(fibs[n][2])) 
  end if

 sequence previous, current, newest --bigfixed
  
    previous = text_num(sprint(fibs[fibslen-1][2]))
    current = text_num(sprint(fibs[fibslen][2]))
    for i = fibslen to n-1 do
            newest = bfi:add(previous, current )
            previous = current
            current = newest
    end for

    return current
end function


--memoize any n>len facts

--**
--eudoc ignore, 
--this really is a private function but want to allow external tests
--ne1 unoptomized version, appoligies to unkmar in advance.
--he provided a few other highly optomized bigfactorials for another project,
--but this is good enough for sporadic usage here. non recursive
--    returns a big num
--this would normally only be used if factorials[] could go over integer size 

export function bigfactorial(integer n) 
  integer lenfacts = length(facts)
  if n <= 0 then
        return text_num(sprint(0)) 
  elsif n <= lenfacts then
        return text_num(facts[n][1]) 
  end if

    sequence b = text_num(facts[lenfacts][1]) --bigfixed 
        --~ n -= 1 

    while n > lenfacts do 
        b = bfi:mul(b , text_num(sprint(n)) )
        n -= 1 
    end while 
  --~ ?123
    return  b  --big num
end function 


function perm_r(object s, atom pnbr, integer len)

  if len = 1 then
    return s
  end if
  
  integer u
    pnbr = remainder(pnbr, factorials[len])
    u = floor(pnbr/factorials[len-1])
    return {s[u+1]} & perm_r(s[1..u] & s[u+2..len], pnbr, len-1)

end function

public function perm_qs(atom pnbr, object s)
  
  return perm_r(s, pnbr-1, length(s))

end function

--**
-- the number of selections of r things from n different things
--  {{{= n! / ( (n-r)! * r! )}}}
public function nCr(atom n, atom r)
  if n - r <= 0 then
      return 1
  --~ elsif n - r < 1 then   --leftover
      --~ return factorials[n]
  end if
  return factorials[n] / ( factorials[n-r] * factorials[r])
end function




--**
--The general formula for the number of combinations 
--of r things chosen from n different things is
--{{{
--
--   .            n!
--   . C(n,r) = ---------
--   .         r!(n-r)!
--}}}
-- from http://mathforum.org/dr.math/
-- not limited  to n<=13 which may be a very large number returned as a string
-- large factorials can be many thousands of digits. 
-- do not inadvertently print results to a terminal you intend to use. (BTDT)
--
public function bignCr(atom n, atom r)
  if n - r <= 0 then
      return text_num(sprint(1))
  end if

  
  --~ return bigfactorial(n) / (bigfactorial(r) * bigfactorial(n-r)
  return bfi:div(bigfactorial(n), bfi:mul(bigfactorial(r), bigfactorial(n-r)) )
end function

--n-combinations of an x-set with repetition allowed.
--x+n-1/n

--~ 9C4 is  --------------- = 126
--~ 35C8 = 32267667
--~ If repetitions are allowed then the calculation is simpler, since we 
--~ can choose every position of a combination of size, say 3, in 35 ways, 
--~ giving 35^3 combinations. For combinations of size 1 to 8, the total is
 --~ 35 + 35^2 + 35^3 + 35^4 + 35^5 + 35^6 + 35^7 + 35^8
 --~ = 2.31810702...  x 10^12
--~ -Doctor Anthony,  The Math Forum
 --~ Check out our web site!  http://mathforum.org/dr.math/  

--isnt hex just a combination problem, 0..F max combinations no repeat $FFFF+1?


--**
-- the number of permutations of r things from n things
--  {{{= n! / (n-r)!}}}
public function nPr(atom n, atom r)
  if n - r < 1 then  --seems like it should be <=1 but that allows some duplcate
      return factorials[n]
  end if
  return factorials[n] - factorials[n-r]
end function


--**
-- the number of permutations of r things from n things
--{{{
--  .                                        n!
--  . The general formula for  P(n,r) is  ------
--  .                                      (n-r)!
--}}}
-- from http://mathforum.org/dr.math/
-- not limited  to n<=13 which may be a very large number returned as a string
public function bignPr(atom n, atom r)
  if n - r < 1 then 
      return bigfactorial(n)
  end if
  --~ return bigfactorial(n) / bigfactorial(n-r)
  --~ return  bfi:div(bigfactorial(r), bigfactorial(n-r)) 
  return  bfi:sub(bigfactorial(r), bigfactorial(n-r)) 
end function


-- The number of partitions of n into k parts
-- not limited  to n<=13 
--~ public function Pk(atom n, atom k)
  --~ if n < k and k = 1 return 0
--~ n (mod k!),

  --~ return 0
--~ end function

--**

-- the number of permutations of combinations of r things from n things
-- not limited  to n<=13 
--unknown if this is accurate yet.
public function bignPCr(atom n, atom r)
  return bfi:mul(bignPr(n, r), bignCr(n, r))
end function



--**
-- the number of derrangements of  n things
--  {{{=  D(n) = n! - n!/1! + n!/2! - n!/3! + n!/4! - ... + (-1)^n*n!/n!}}}
-- http://mathforum.org/dr.math/
--the number of derangements can be approximated by n!/e
--{{{
--!N = No.of ways of arranging N things such that none of them appear in their natural order.
--The general formula for this is:
--!N = Sum(r=0 to n) [(-1)^r * nCr * (n-r)! ]
--}}}
--
public function nD(atom n)
  --~ sequence D={0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961}
  sequence D={
     1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961
     , 14684570, 176214841, 2290792932, 32071101049
     --~ , 481066515734, 7697064251745, 130850092279664
     --~ , 2355301661033953, 44750731559645106
     --~ , 895014631192902121, 18795307255050944540  
   }
  if n < 1 or n>length(D) then 
    printf(1, "nD(n)  unfinished for n=%d\n", {n} )
      return 0
  end if
  return D[n]
end function


--~ Unsigned Stirling numbers of the first kind
    --~ c(n,k)=\left[{n \atop k}\right]=|s(n,k)|=(-1)^{n-k} s(n,k)
--~ (with a lower-case "s") count the number of permutations of n elements with k disjoint cycles.
--~ The Stirling numbers of the first and second kind can be considered to be inverses of one another:
--~ function S(atom n, atom k)
  --~ sequence S={}
    --~ return S[n][k]
--~ end function



-- S(n,k) Stirling numbers of the second kind 
--~ function S(atom n, atom k)
  --~ sequence S={}
   --~ return S[n][k]
--~ end function
 

--**
-- The Bell number B(n)  is defined as 
-- the number of partitions (of any size) of a set 
-- into n distinguishable objects.
-- http://www.research.att.com/~~njas/sequences/A000110
-- http://en.wikipedia.org/wiki/Bell_numbers
public function Bn(atom n)
  sequence Bells={
        1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147
        , 115975, 678570, 4213597
  , 27644437 , 190899322, 1382958545 , 10480142147
  --~ , 82864869804, 682076806159, 5832742205057
  --~ , 51724158235372, 474869816156751
  --~ , 4506715738447323
   }

  if n < 0 or n-1>length(Bells) then 
      printf(1, "Belln(n)  unfinished for n=%d\n", {n} )
      return 0
  end if
  
  return Bells[n-1]
end function

--~ -- 
--~ public function parity(atom n)
  --~ return sum(rs)%2
--~ end function


function slct(object s, atom n, atom r, atom q)

  -- this recursive function does the real work

  integer h

  if r = n then
    return s
  elsif r = 1 then
    return {s[q]}
  end if

  h = nCr(n-1, r-1)

  if q <= h then  -- 1st element is s[1], select r-1 from s[2..n]
    return {s[1]} & slct( s[2..n], n-1, r-1, q)
  end if
     -- select r from s[2..n]
    return slct( s[2..n], n-1, r, q-h)

end function

--**
-- the qth selection of r from s
public function select_qrs(atom q, atom r, object s)

  -- this calls slct() which does the real work

  return slct(s, length(s), r, 1+remainder(q-1,nCr(length(s),r)))

end function

function indices(atom n, atom r)
  return {1+remainder(floor((n-1)/r),r), 1+remainder(n-1,r)}
end function

--**
-- calculate the qth permutation of r elements from s
public function perm_qrs(atom q, atom r, object s)
    
  object s2, ix

  ix = indices(q, factorials[r])
  s2 =  select_qrs(ix[1], r, s)
  return perm_qs(ix[2], s2)

end function



--**
-- added as a convenience
-- return all combinations
-- of r elements from s,  make tosort FALSE to return as is
-- defaults to RD_SORT return sequence sorted with only unique items.
-- set to RD_INPLACE just remove duplicates. ##include std/sequence.e##
-- 
-- r defaults to length of s

public function combs(sequence s, integer r=0, integer tosort=RD_SORT)
  sequence results={}, t
   if r <= 0 then r = length(s) end if

    for j = 1 to nCr(length(s), r) do
	    results = append(results, select_qrs(j, r, s))
    end for

  if tosort then  --
    return remove_dups(results, tosort)
  end if
  
  return results
end function



--**
-- added as a convenience
-- return all the permutations of all the select combinations
-- of r elements from s,  make tosort FALSE to return as is
-- defaults to RD_SORT return sequence sorted with only unique items.
-- set to RD_INPLACE just remove duplicates. ##include std/sequence.e##
-- 
-- r defaults to length of s
--
-- Caution: for more than 4 characters this can be a very large sequence
-- there probably is loss of efficiency calling for the nth permutation
--instead of looping inline directly, that may happen in a future version.
-- paste this or similar code in your program to avoid handling large sequences.

public function permcombs(sequence s, integer r=0, integer tosort=RD_SORT)
  sequence results={}, t
   if r <= 0 then r = length(s) end if

    for j = 1 to nCr(length(s), r) do
          t = select_qrs(j, r, s)
	for k = 1 to nPr(length(t), r) do --1..length(t) not enough
	    --~ results &= {perm_qs(k, t)}
	    results = append(results, perm_qs(k, t))
	end for
    end for

  if tosort then  --
    return remove_dups(results, tosort)
  end if
  
  return results
end function


--**
-- same as select_qrs routine but limited to r under 5 to avoid long startup time
--uses precalculated sequence of all permutations and combinations indexes
public function select_qrsm(atom q, atom r, object s)
  return 0
end function

--**
-- same as perm_qs routine but limited to r under 5 to avoid long startup time
--uses precalculated sequence of all permutations and combinations indexes
public function perm_qsm(atom q, object s)
  return 0
end function


--**
-- same as permcombs routine but limited to r under 5 to avoid long startup time
--uses precalculated sequence of all permutations and combinations indexes
--then translates your sequence using those indexes to save time
--on a small sequence this probably will not be noticable
--but will allow testing on larger sequences once tested on smaller sequences
public function permcombm(sequence s, integer r=0, integer tosort=RD_SORT, integer uniques=0)
   --need to get a mapping of nCr nPr to see how many precalculated will need
   --~ if r > 5 or length(s) > 10 then 
   if nCr(length(s), r) > 6 or length(s) > 5 then 
       printf(2, "your r or s is too large for this routine\n",{})
       return s 
   end if
   
   sequence results={}, t
   if r <= 0 then r = length(s) end if

    for j = 1 to nCr(length(s), r) do
          t = select_qrsm(j, r, s)
	for k = 1 to nPr(length(t), r) do --1..length(t) not enough
	    --~ results &= {perm_qs(k, t)}
	    results = append(results, perm_qsm(k, t))
	end for
    end for

  if tosort then  --
    return remove_dups(results, tosort)
  end if
  
  return results
end function



with trace
--~ trace(3)

--**
--emulate python destructive pop from end of array, 
--returns the sequence.
--for the converted pynms functions 

export function pop_seq(object seq)
   if not sequence(seq)  then
	  ? 100 & seq
		return ""
  elsif  length(seq) =1 then
	 if sequence(seq[1]) then
		return {}
	 else
		return ""
	 end if
   elsif not length(seq) then
	  ?3333
		return {}
   end if
  return seq[1..$-1]
end function

without trace



--**
--emulate python pop from end of array, 
--call pop_seq() to remove last item.
--for the converted pynms functions 
--
--seems to be working either string or sequence of atoms or sequences
-- problematic, is 'z' the last element of "xyz" still a sequence 
--.still wrapping in braces sometimes.need more exaustive tests 

export function pop_tail(object seq)
   if not sequence(seq)  then
		return {seq}
   --~ elsif  length(seq) =1 then
	 --~ if sequence(seq[1]) then
		--~ return seq
	 --~ else
		--~ return {seq}
	 --~ end if
   --~ elsif not length(seq) then
	  --~ ?-999
		--~ return {}
   end if
  return seq[$]
end function


--**
-- convert from pynms 
--none functional yet
--
--
-- convert from pynms 
-- algorithm by Christos Georgiou 
-- Returns the <index>th permutation of <seq>
--
--converting get_* from python a little tricky, objects are not by reference
--so changes inside the function are not reflected to each recursion.
--transparently handling difference between {1} and "1" not simple
--still producing too many duplicates and not covering every posibility

public function get_permutation(object seqc, atom index)
    --~ object seqc = seq   -- list(seq[:])
    object seqn =  pop_tail(seqc)    -- seqc.pop()
     seqc = pop_seq(seqc)    -- finish pop()
    integer divider = 2 -- divider is meant to be len(seqn)+1, just a bit faster
    integer new_index
    while sequence(seqc)  and length(seqc) do
        --~ index, new_index= index//divider, index%divider
        index = floor(index/divider)
        new_index = remainder(index, divider)
	  --~ ?seqn
	    if sequence(seqn) and length(seqn) <= new_index then 
			seqn = insert(seqn, pop_tail(seqc), new_index)   -- seqc.pop()
		else
		  seqn &= {pop_tail(seqc)}
		end if
        seqc = pop_seq(seqc)  --finish pop
        divider += 1
    end while
  return seqn

end function

with trace
--~ trace(3)

--**
-- convert from pynms 
-- # XXX
--
-- currently both combination routines are returning the empty set
--
public function get_unique_combinations(object alist, atom numb, object blist={})
	if not numb then
		return blist
    end if
	sequence o
	--~ for i in range(len(alist)) do
	for i = 1 to length(alist)-1 do
            blist = append(blist, alist[i])
	        --python may be handling default arg and blist differently
		o = get_unique_combinations(alist[i+1..$], numb-1, blist)
		alist = o[1]  blist = o[2] 
        blist = pop_seq(blist)  -- blist.pop()

	end for

	return {alist, blist}   --> a="a" b="b" o="ab" woops

end function

without trace

--**
-- convert from pynms 
--
-- currently both combination routines are returning the empty set
--
public function get_combinations(object alist, atom numb, object blist={})
	if not numb then
		return {alist, blist}
    end if
	sequence o
	--should you iterate over a list whiile changing it?
	for i = 1 to length(alist) do
		blist = append(blist, alist[i])  -- alist.pop(i))
	    --probably should to range check
        alist = remove(alist, i) --finish pop
		o = get_combinations(alist, numb-1, blist)
		alist = o[1]  blist = o[2] 
		alist = insert(alist, pop_tail(blist), i)  -- blist.pop()
        blist = pop_seq(blist) --finish pop
    end for
  return {alist, blist}

end function

--r09827p05:31
--eucode is barfing on something badsig at perm_r and skips the end routines
--doesn't like the formula for nCr


--get_perms 
--have to test it to make sure it can return every posible permutation
--seems to wrap on q past possible and return many duplicates
--
-- both get combination routines are returning the empty set
--uncle
--
--renamed nCr() to  nPr and added nCr()
--simplified some code
--not return factorials[n] when n-r=1
--
-- fixed, why is it named select_rs if it's really qs?
-- why not call them comb_? instead of select_?
--need to allow repitions with an option off
--not to be confused with duplicates which may be a bug
--still testing. hard to find authoritive exaustive non trivial test cases
--could still be mixing up char string and numeric sequences
--possibly subsequences will not combine well or even too well
--turning a simple permutation into matrix math
--added permcombs routine to get all r at a time from string
-- fixed nCr to return 1 when n-r =0 instead of factorials[n]
-- seems like only one selection is possible when r=n
--posted permcombs nPr and fixed nCr to forum
--chg results &= to append in permcombs()
--derek fixed a leak that could have impacted the test on very large permutations
-- results[n] &= something. I didn't notice a problem with perms6
--and it was reassigning after each test. should run again and check usage
--

--possibly should be included in the stdlib with a few more changes and bug fixes?
--maybe drag in bigint for longer than 13!

--add a compressed format for returns from permcombs
--helper functions to locate an index and to/from the compressed mode
--could even be simple runlength encoding
--eventually to/from a db with code to regenerate the db to a max size or beyond

--for combinations and permutations less than 5, can actually include directly
--would need a little cleaner way to access them than used in the test program
--could even memoize, euphoria needs a gp memoize routine
--maybe there is one in primes already?
-- could use a decent timeit program, maybe there is on in eu4/bin now?
--run a unittest with larger values and see if it shows up in the html report


--maybe nCr nPr is not optomal either
--there was a forum post about duplicate letters with no repitition
--returning many less permutations than one where every item is unique
--maybe should at least compare a uniqued string to a permuted then uniqued
--to see if can pre unique it and have the same results

--why did I switch to object? is there any speed peanelty?
--need to compre sequence & object interpreted & translated
--with watcom and minGW in some kind of longer then a minute test

--add a random filled array, expecially for testing
--specify the size minimal distance or somthing, 
--and it skips so many and either zero fill or condense
--this could allow much larger permutations than time would otherwise allow

--added functions for factorial >13
--~ Factorials with prime factorization

--add append if > len of existing sequence
--not sure that is threadsafe used in a dll though


